Microscopic Current Dynamics in Nanoscale Junctions 
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So far transport properties of nanoscale contacts have been mostly studied within the static scattering ap- 
proach. The electron dynamics and the transient behavior of current flow, however, remain poorly understood. 
We present a numerical study of microscopic current flow dynamics in nanoscale quantum point contacts. We 
employ an approach that combines a microcanonical picture of transport with time-dependent density-functional 
theory. We carry out atomic and jellium model calculations to show that the time evolution of the current flow 
exhibits several noteworthy features, such as nonlaminarity and edge flow. We attribute these features to the 
interaction of the electron fluid with the ionic lattice, to the existence of pressure gradients in the fluid, and to 
the transient dynamical formation of surface charges at the nanocontact-electrode interfaces. Our results sug- 
gest that quantum transport systems exhibit hydrodynamical characteristics which resemble those of a classical 
liquid. 



I. INTRODUCTION 

Recent experimental progress has enabled imaging of co- 
herent current flow dynamics in quantum point contacts 
formed in semiconductor heterostructures i *' 2 ' 3 ' 4 These ad- 
vances in experimental techniques open the possibility that 
current flow through atomic or molecular junctions will be 
eventually imaged and controlled. Understanding the micro- 
scopic electronic flow patterns can aid the design of novel 
electronic devices. However, very few theoretical studies of 
current dynamics in nanoscale systems are currently avail- 
able. Indeed, among the recent theoretical studies of trans- 
port in nanoscale systems, much emphasis has been placed 
on the steady-state conduction properties ^ 5 ' 6 ' 7 ' 8 ' 9, 10 ' 'WW 5 
whereas the transient behavior of the current remains an un- 
explored area. 

Electronic transport in nanoscale junctions is usually for- 
mulated within the stationary scattering picture, such as the 
one due to Landauer, 16 in which the conduction is treated as 
a collection of scattering events. This stationary approach, 
widely used to study transport in mesoscopic and nanoscale 
systems, has led to considerable success in understanding 
current-related effects other than the conductance, including, 
e.g., noise, local heating, and current-induced forces. 17 It has 
helped our understanding of the microscopic current distri- 
bution as well4&i2 Nevertheless, the stationary approach as- 
sumes that the system is already in a steady state, leaving the 
questions of how a steady current establishes itself and what 
other phenomena are related to the dynamical formation of 
steady states in a nanojunction unsolved. 

The dynamical nature of the current flow is better addressed 
in a time-dependent approach than in the stationary one. 
Time-dependent or AC transport approaches have been pre- 
viously introduced in mesoscopic conducting systems . 2Q i 21 i 22 
Recently, an increasing amount of effort has been directed 
toward developing ab-initio time- dependent approaches for 
nanoscale systems . 23 ' 24 ' 25 ' 26 ' 27 ' 28,29 ' 30 ' 31 ' 32 ' 33 One such method 
was developed to study the AC conductance using the time- 
dependent density functional theory 34 (TDDFT) combined 



with absorbing boundary conditions. 24 This method, however, 
is affected by the arbitrariness of how and where the ab- 
sorbing potentials are added, while the effect of the absorb- 
ing potential on the conduction in nanojunctions is unclear. 
Other methods have been developed based on the Landauer 
scattering formalism of transport that employ open boundary 
conditions, 25 ' 27 ' 29 ' 30 

More recently, a microcanonical formalism that treats elec- 
tronic transport as a discharge across a nanocontact connect- 
ing two large but finite charged electrodes was introduced. 23 
The formulation has been combined with TDDFT to study the 
dynamical formation of quasi-steady currents, local chemical 
potentials, 28 and electron-ion interactions. 26 This formalism 
would yield the exact total current flowing from one electrode 
to the other if the exact functional were known, regardless of 
whether the system achieves a steady state. In practice, in ab- 
initio transport calculations, one only uses approximate forms 
of the functional such as the adiabatic local density approx- 
imation (ALDA). It has been shown recently that the elec- 
tronic correlation effect beyond the ALDA gives rise to ad- 
ditional resistance in molecular junctions i 35 ' 36 ' 37 The spurious 
self-interaction implicit in the ALDA further complicates cal- 
culations of the conduction properties. 39 ' 40,41 The sensitivity 
of various transport properties to the suggested corrections re- 
mains only partially understood. 

The dynamical establishment of a quasi-steady current has 
been investigated by a number of authors i 26,27,28 ' 42 It has been 
demonstrated that a quasi-steady current can establish itself 
across a junction on a femtosecond time scale without the 
presence of inelastic scattering. 28 This is due to the geomet- 
rical "squeezing" experienced by the electrons crossing the 
nanojunction. 23 The conductance calculated from the micro- 
canonical formula was shown to be in good agreement with 
that obtained from the static DFT approach in prototypical 
atomic junctions 26,28 as well as in molecular junctions. 43 Nev- 
ertheless, a study of the microscopic behavior of the electron 
flow, and in particular of the current flow morphology in nano- 
junctions, is still lacking. 

In this paper, we carry out real-time numerical simulations 
of current flow in metallic nanojunctions using the micro- 



2 



X 



V(z); 



Za 



FIG. 1: (Color online) Sketch of the nanojunction geometry that 
is studied in the paper. At t < 0, a bias in the form of V(z) — 
Vb[H(z — z a ) — H(— z — z a )] is applied to the junction (the central 
constriction is at z — 0) such that the regions \z\ > z a bear a poten- 
tial offset from the central constriction, where H(z) is the Heaviside 
step function. 



canonical formalism, where we employ TDDFT within the 
ALDA. We restrict the forthcoming discussion to the dynami- 
cal behavior of electron/hole charges in the nanojunctions un- 
der the linear response regime, i.e., in which the bias is small 
and the current-voltage characteristics are linear. The paper 
is organized as follows. In Sec|]]]we discuss model transport 
systems and numerical methods. In SecJTTTJ we present and 
discuss simulations of current dynamics in jellium and atomic 
junctions. We also analyze the effects of hydrodynamic pres- 
sure and electrode surface charges on the dynamics of the 
flow. In SeclIVI we summarize the main conclusions of our 
work. 



II. MODEL AND METHODS 

The nanoscale junction geometry studied in this paper is il- 
lustrated in Fig.Q] A narrow constriction separates two large 
but finite electrodes. We begin the simulations by applying a 
step function-like electric bias across the junction such that the 
two electrodes bear equal and opposite potentials offset rela- 
tive to the potential at the center of the junction. The distance 
from the bias discontinuity to the center of the junction is z a 
(see FigQ]). This bias induces a charge imbalance between the 
two sides of the system. At t = 0, we remove the bias and a 
discharge through the nanojunction ensues. The Kohn-Sham 
initial state therefore corresponds to the ground state of the 
Hamiltonian in the presence of the bias. Here, we are inter- 
ested in the transient behavior during the phase in which the 
current is in the process of establishing a quasi-steady state 
and immediately thereafter, i.e., long before the electrons that 
have passed through the constriction have had chance to reach 
the far boundary of the electrodes. 

To separate the effects of the electrons and of the atomic 
lattice, we have carried out calculations using a jellium model 
and an atomic model. In the jellium calculations, the elec- 
trodes are represented by two large jellium slabs 2.8 A thick. 
The contact is a rectangular jellium block 2.8 A wide and as 



thick as the electrodes. The distance between the jellium elec- 
trodes is 9.8 A. In the atomic calculations, the junction is rep- 
resented by two planar arrays of gold atoms sandwiching a 
single gold atom. We employ TDDFT and solve the effective 
single-particle Schr odinger equation 



ihip(f, t) = 



h 2 v 2 

2m 



+ V c g(r,t) 



(1) 



where the effective potential is given by 



V cB (r, t) = V ext (f ; t) + / ^^df + 14c (r , t) . (2) 



The term V xc (f, t) is the exchange-correlation potential cal- 
culated within the adiabatic local- density approximation. The 
external (ionic) potential is modeled using pseudopotentials 
for the atomic calculations, 44 while in the jellium model it is a 
local operator related to the uniform positive background jel- 
lium density po via V cyi t(r) = J -J^dr 1 , where po equals 

to inside the jellium and zero outside, and r s is 

the Wigner-Seitz radius. In the jellium model, we choose 
r s = 3a,B (as = Bohr radius) which gives a good repre- 
sentation of bulk gold (see also discussion below). A "free- 
space" boundary condition is implemented such that the long- 
range potential is constructed only from the densities in the 
supercell; 45 that is, the system is not periodic. Additional 
numerical details can be found in Appendix [A] The single- 
particle time-dependent current density is calculated via 
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r n (f, t)V^„(r, t) - V<(f, t)Mr, t) 



(3) 

where ip n denotes individual Kohn-Sham single-particle 
states. Note that in TDDFT the current density is not neces- 
sarily exact even if one calculates it with the exact functional 
(via the continuity equation, only the divergence of the current 
density would be exact). One would need to resort to Time- 
Dependent Current Density Functional Theory (TDCDFT) 46 
to obtain the exact current density (with the exact functional). 
Nonetheless, due to the small viscosity of the electron liquid, 
we have found that the current density one obtains by using 
TDDFT within the ALDA, and the one obtained by TDCDFT 
within the Vignale-Kohn functional 46 are qualitatively simi- 
lar41 

Even if the contact between the electrodes were removed, 
the current between the two electrodes would not completely 
vanish because of quantum tunneling. This bare tunneling 
current can conveniently be used to compare the jellium and 
the atomic calculations. The jellium edges are placed at half 
the interplanar spacing of the lattice. 48 This way, the jellium 
model and the pseudopotential calculations both yield tunnel- 
ing current densities of ~ 0.05 /1A/A 2 at a bias of 0.2V. The 
agreement indicates that the jellium model is a good represen- 
tation of two large metal electrodes. This is consistent with 
the results of previous density-functional calculations. 49 
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FIG. 2: (Color online) Current flux for a series of times in a nanoscale quantum point contact system in the jellium model. The applied bias at 
t < is A V — 0.2V. The field lines in each panel depict the direction and amplitude of the current density vectors, while the colors give extra 
electron (red) or hole (blue) density, (a) t = 0.4 fs ; (b) t = 0.8 fs ; (c) t = 1.6 fs. In (a), the semicircle marks the contour along which the 
radial component of the current density is calculated (see text and Fig.[5J. 



III. RESULTS AND DISCUSSION 
A. Flow dynamics through jellium model junctions 

In a nanojunction such as an atomic point contact, the di- 
mensions of the leads are usually much larger than those of 
the central constriction. In addition, not far from the con- 
tact, we expect the electron momentum to converge to the 
value characteristic of the bulk leads. Therefore, the momen- 
tum of an electron coming from the leads and entering the 
contact has to change considerably. This gives rise to resis- 
tance, and for a truly nanoscale junction, this momentum mis- 
match is mainly responsible for the establishment of quasi- 
steady states i 23 ! 28 ' 50 Using the above dynamical approach we 
can now study how this translates into microscopic current 
flow through the nanocontact and into the leads by calculating 
the current density at different times. 

To begin, we follow the method described in the preced- 
ing section to impose a charge imbalance in the jellium model 
system. A discussion of the effect of the lattice on the flow dy- 
namics will be presented in the following section. The initial 
bias is chosen such that the discontinuity happens at the edge 
of the jellium slab near the central constriction (i.e., z a = 5 
A). The flow pattern is independent of the location of the dis- 
continuity once the current starts to flow through the center of 
the junction. In Fig. [2] we plot three snapshots of the current 
density to illustrate the evolution of the flow. Due to the bias 
offset near the jellium edges, a dipolar layer forms on each 
of the two contact-electrode interfaces. As a result, the initial 
current flow is uniform on both sides as shown in Fig. EJa). 
Very little current flows in the nanojunction at this point, how- 
ever the current steadily rises. In Fig.[2fb), the current density 
becomes convergent toward the center of the nanojunction. In- 
terestingly, as the excess charge from the left electrode reaches 
the contact, there is a period of adjustment during which the 
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FIG. 3: (Color online) Time series of the radial amplitude of the cur- 
rent densities along a semicircle of radius 3A centered on the nano- 
junction, as a function of angle along the contour. 



dominant flow is in the lateral direction, i.e., parallel to the 
facing surfaces of the electrodes. 

To quantify the evolution of the angular distribution of the 
electron flow, in Fig.[3]we plot a time series of the radial com- 
ponent of the current density along a semicircle contour cen- 
tered on the junction as a function of the angle on the semi- 
circle (see Fig. |2] (a)). One can see that initially a "wave" of 
excess charge approaches the nanocontact. Then, the radial 
current density peaks at very large angles (~ ±75°); i.e., the 
current density near the contact is dominated by the flow along 
the electrode edges in the lateral direction. The peaks then 
gradually move towards the central axis, and the current den- 
sity adjusts to a more "focused" pattern as shown in Fig. [2c). 
We have also examined a junction consisting of a jellium cir- 
cular "island" between the electrodes. We have also observed 
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FIG. 4: (Color online) Time sequence of electron current streamlines in the atomic junction described in Sec|n] Panels (a) - (d) correspond 
to t = 0.2, 0.25, 0.3, and 0.35 fs, respectively. The dots denote atomic sites that corresponds to the (001) facets of the gold FCC lattice. The 
applied bias at t < is AV = 0.2V. 



edge flow in this case as well. The edge flow is not a quantum 
interference pattern, and cannot be compared with the fringes 
observed in the 2D electron gas quantum point contacts^ In- 
stead, we suggest that the flow pattern is controlled by hydro- 
dynamic effects and forces due to surface charges. We analyze 
these in Sec lHICI 

The structure of the current density is analogous to a clas- 
sical fluid flowing across a narrow constriction. This is not 
surprising because an inhomogeneous electron system can 
be indeed characterized by a set of hydrodynamical rela- 
tions expressed in terms of the particle density and velocity 
field ! 51 ' 52 ' 53 ' 54 More recently, in particular, a hydrodynamical 
approach was proposed for nanoscale transport systems^ fur- 
ther strengthening the analogy between the dynamics of the 
electron liquid and the one of a classical liquid. We note that 
the present calculations do not take into account the physi- 
cal viscosity of the electron liquid. 54,55 An inviscid fluid can 
therefore be used as a model for the dynamical behavior of the 
present electron liquid. 



B. Flow dynamics through atomic junctions 

The jellium model was a convenient way to probe the mi- 
croscopic current dynamics in an electron gas. To understand 



the influence of the lattice on the flow, we carried out simu- 
lations that included an ionic background modeled by pseu- 
dopotentials. The atomic calculations were carried out for 2- 
D gold nanojunctions and were initialized in the same way 
as the jellium calculations. The theoretical and experimental 
conductance are in very good agreement for this system^ We 
chose lattice arrangements that correspond to the (001) and 
(111) facets of the gold FCC lattice (the dots in Fig. |4]repre- 
sent the atomic sites for the (001) configuration). 

Electric current streamlines at different times in the simu- 
lation are plotted in Fig. HJa) - (d). The streamlines are cal- 
culated by integrating the current density field upstream and 
downstream, dr/ds = ±j(r(s)). The morphology of the cur- 
rent flow in the atomic junction and in the jellium model is re- 
markably similar, 57 indicating that the jellium model is a good 
representation of the gold electrodes. Nevertheless, a number 
of new features appear in the atomic calculations. 

Fig-Hfc) shows that once a steady flow through the junction 
is established, the current spreads into a wedge-shaped region 
inside the electrodes. The flow morphology for each of the 
two different lattice arrangements is similar except that the 
flow spreads over a broader wedge-shaped region in the (111) 
lattice. Another common feature in the atomic calculations is 
the presence of a stagnant zone at the corner of the electrode 
boundary. There is little current flow into or out of this zone. 



This is similar to a classical fluid where a stagnant zone can 
sometimes be located at the entrance or exit of a channel. 

One profound difference between the atomic and the jel- 
lium calculations is the formation of eddies evident in the for- 
mer but not in the latter. In the jellium calculations carried out 
within the linear-response bias regime, the current flux lines 
are laminar. In contrast, in the atomic calculations, the ed- 
dies appear as localized circular flow that can be observed in 
Fig. Hfd). The eddies develop in both electrodes and the size 
of the eddies is comparable to the interatomic distance. The 
eddies are reminiscent of the vortices that form in a classical 
fluid at higher Reynolds number when the fluid encounters 
obstacles. As is well known, vortices can also occur when 
velocity shear is present within a continuous fluid (Kelvin- 
Helmholtz instability). We suggest that the lattice ionic ob- 
stacles and the boundaries separating the flow zone and the 
stagnant zone facilitate the formation of the observed eddies 
in our simulations. 

The formation of current vortices has been previously re- 
ported in 2D ballistic quantum billiards. 58 ' 59 In these quantum 
systems, a rich variety of flow patterns ranging from regular 
to chaotic is possible. While we cannot draw direct analo- 
gies with these open and mesoscopic transport systems, an 
unstable and turbulent flow has been recently predicted in 
nanoscale systems 53 . We have argued that the ALDA elec- 
tron liquid in our simulations corresponds approximately to 
an inviscid fluid. We therefore suggest that in the presence of 
a lattice, hydrodynamical instabilities, or turbulence can oc- 
cur in nanotransport systems even in the absence of a physical 
viscosity. 



C. Hydrodynamics and the formation of surface charges 

To understand what drives the edge flow along the electrode 
surfaces, we examine the evolution of the charge distribution 
near the surfaces of the junction. For this purpose, we ap- 
ply a step-function bias such that the potential discontinuity 
in each electrode occurs at z a = 10A, which cuts across the 
electrode, as illustrated in Fig.Q] In Fig.|5ja), we plot a time 
series of the x — y plane-averaged excess charge density along 
the z axis. At t = 0, two symmetric dipolar layers form in- 
side the electrodes as a result of the bias offsets. As the cur- 
rent starts to flow through the contact and gradually reaches a 
quasi-steady state, a global charge redistribution becomes ap- 
parent. The dipolar layers gradually vanish and are replaced 
by surface charge layers that form at the contact- electrode in- 
terfaces. The charge contour plots Fig.|5jb-c) further illustrate 
the formation of surface charges as a result of current flow. 

The formation of surface charges around the central con- 
striction is reminiscent of the formation of residual-resistivity 
dipoles introduced by Landauer. 16 It has been suggested that a 
continuous current flow arriving at a junction must be accom- 
panied by self-consistently formed charges at the electrode 
surfaces. 49 The effect should be taken into account to correctly 
characterize the electrostatic potential and the nonequilibrium 
conducting properties in a transport calculation. 61 In this work 
we provide the first numerical demonstration of the dynami- 
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FIG. 5: (a) Planar averaged charge density from t = to t = 0.6 fs 
for a jellium junction. The change in the peaks indicates the dynamic 
process in which excess charge builds up at the surfaces. The sign 
of the surface charges indicates that electron charges accumulates on 
the right and hole charges on the left, (b) - (c) Excess charge in the 
vicinity of the contact at t = and t = 0.6 fs. Thick solid lines 
indicate excess electrons, while thick dashed lines indicate excess 
holes. Thin dashed lines mark the edges of the jellium electrodes 
and the contact. 



cal formation of these surface charges using a time-dependent 
approach. 

We have already observed that, as the surfaces of the elec- 
trodes are populated by excess charges, a lateral flow starts 
to develop along the surfaces. This behavior is illustrated in 
Fig. [3] where the radial current flux at t = 0.8 fs shows two 
pronounced peaks at very large angles. We attempt to inter- 
pret this behavior within the framework of an effective clas- 
sic hydrodynamic model of an inviscid charged fluid. The 
acceleration of the fluid is then given by Euler's equation 
dv/dt + (v ■ V)v = — VP/m e n — eV<p/m e , where v is the 
fluid velocity, n is the fluid particle density, m e is the elec- 
tron mass, and ip is the electrostatic potential. The first term 
on the right hand side is the acceleration due to a gradient in 
electron pressure. The second term is the acceleration due to 
the electric field of the excess charges on the surfaces of the 
electrode. The electric field drives the electrons/holes toward 
the surfaces to cancel out the excess charges. 

The inertial term in the above equation can be estimated 
as \(v ■ V)v\ ~ v 2 /L ~ 10 6 m 2 /s 2 x L~ x . Here, L de- 
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notes the characteristic length scale over which we expect a 
departure from uniform flow, so that V ~ L^ 1 . Velocity 
of the flow in the simulation reaches v ~ 10 3 m/s <C vjr, 
where Vp ~ 1.4 x 10 6 m/s is the Fermi velocity of bulk 
gold. The hydrodynamical pressure can be calculated from 
the derivative of the ground state energy with respective to r s , 
P/n — — ^e'(r s ), where e(r s ) is the energy per particle. 62 

Therefore, \VP\/m e n ~ ^v F \Vn\/n. Here, we included 
only the ground state energy of a noninteracting electron gas, 
£o( r s) = §£f-— Let Sn denote the change of the particle den- 
sity as a result of the current or the formation of the surface 
charges. For the present junctions at these biases, we find that 
typically Sn/n < 0.01. Then, the acceleration due to the pres- 
sure gradient is 

_ |VP| lv 2 F Sn 

\ap\ = ~ - — 

m e n 5 L n 

< 3 x 10 9 ^ x L-\ (4) 

To estimate the magnitude of the electrostatic acceleration, 
we treat the layer of charge induced on the facing surfaces 
of the electrodes (as illustrated in Fig[5J;) as an infinite uni- 
formly charged wire. The electric field of the wire is given 
by |Vy| = \I1-keqL, where A is the linear density of excess 
charge, L is the distance to the wire, and £o is the vacuum per- 
mittivity. The linear charge density is calculated by averaging 
the charge density difference e8n between the configuration 
with and without the current flow over the layer in which the 
charge accumulated at the contact-electrode interfaces. For 
the same junctions, we find A ~ 0.016 e/A. The acceleration 
of charges due to this electric field can be calculated as 



„ 5xl0 11 ^xr 1 . 
These crude estimates imply that over the same distance L 

\(v- V)v\ < \a P \ < \S d \, (6) 

which suggests that the hydrodynamic pressure gradients 
dominate over the inertia of the fluid (the flow is subsonic and 
compressible), while the maximum electrostatic force due to 
the surface charges is comparable to or larger than the pres- 
sure gradient force before the surface charge has been passi- 
vated. Therefore, it is plausible that the lateral flow observed 
in the simulation is primarily of electrostatic origin. In differ- 
ent junction geometries or in a different conductance regime, 
the ordering in Eq.|6]may be different. 

We have also carried out a similar simulation using a 
parabolically- shaped constriction that resembles a quantum 
point contact in the 2D electron gas. At a similar bias as in the 
non-parabolic junctions, we find similar surface charge accu- 
mulation along the boundaries of the electrodes in the vicinity 
of the contact. We believe the analysis we have provided on 
the surface charge formation applies to this case as well. To 
the best of our knowledge, the accumulation of the surface 



charges has not been reported in adiabatic quantum point con- 
tacts before. It would thus be interesting to develop experi- 
mental techniques to explore the surface region of the quan- 
tum point contact in a 2D electron gas and the charge accu- 
mulation that we observe in our simulations. 



IV. CONCLUSIONS 

In this paper, we have used the microcanonical approach 23 
to study the time-dependent current flow morphology and the 
charge distribution in discharging nanojunctions represented 
by both a jellium model and pseudopotentials. We have 
showed that the electron flow in the nanojunctions exhibits hy- 
drodynamic features analogous to a classical fluid. We have 
found that in the atomic case the current flow evolves into 
wedge-shaped pattern flanked by stagnant zones. The flow 
develops nonlaminar features including eddy currents. We 
suggest that the ionic lattice at the junction plays the role of 
"obstacles" in the dynamics of the electron liquid, with con- 
sequent development of these features. We have also demon- 
strated that excess surface charges accumulate dynamically 
along the electrodes. In addition, we have observed that for 
a period of time, there is strong current flow in the transverse 
direction. By employing an order-of-magnitude argument we 
suggest that this flow is driven by both hydrodynamic forces 
due to the electron pressure, and electrostatic forces due to 
the surface charge distributions. The latter forces dominate 
the initial dynamics in the junctions at hand. 

The present and a previous study 28 demonstrate that 
the microcanonical approach combined with time-dependent 
density-functional theory can be used to probe the transient 
behavior of the current in nanojunctions such as atomic-scale 
point contacts or molecular junctions. 43 The present approach 
supplements existing methods that are based on the static scat- 
tering picture and provides another tool to studying relatively 
unexplored nanoscale transport phenomena from first princi- 
ples. 

The flow patterns we observe in metallic nanojunctions can 
be generalized to a number of other systems, such as molecu- 
lar junctions, although many details may vary. In view of the 
recent advances in microscopic imaging techniques of coher- 
ent current flow in quantum point contacts in a 2D electron 
gasr- we hope that new experimental work exploring the be- 
havior of current flow in atomic contacts and molecular junc- 
tions will soon emerge to supplement our studies. 
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APPENDIX A: NUMERICAL DETAILS 

We performed time-dependent density-functional calcula- 
tions using the program socorrc 44 and an in-house pro- 
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gram which implements TDDFT within the jellium model. 
The gold ions were modeled by norm conserving Hamann 
pseudopotentials with 6s electrons as valence electrons.— We 
used the Perdew-Zunger (1981) LDA exchange-correlation 
functional. 66 The atomic calculation employed a plane-wave 
basis set, with an energy cutoff of 204 eV, which corresponds 
to grid spacing of 0.2A. The energy eigenvalues vary by less 
than 1% by increasing the cutoff by 66%. In the jellium case, 



the calculations were performed using a real-space basis set 
where the space is uniformly discretized and the grid spac- 
ing is 0.7 A. The eigenvalues vary less than 3% by decreasing 
the grid spacing by 66%. The time evolution operator is rep- 
resented using the Chebyshev method, 67 with a time step of 
5 x 10~ 4 fs in the atomic case and 5 x 10~ 3 fs in the jellium 
case. 
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